Data: annotate post-2020 reactions with pathways; port the step that was lost in 2020 - #300

A data defect found while building the pathway panel reviewer 1 asked for (PR #299). Independent of the manuscript — no latex/ changes here.

Withdrawn 2026-09-28 — see the last section. The branch is back on its standalone dev-based commits.

The defect

Not one of the 12,261 live reactions added since the 2020 release (rxn48576–rxn60859) carries a pathway or EC annotation.

Unique_ModelSEED_Reaction_Pathways.txt stops at rxn48568 and Unique_ModelSEED_Reaction_ECs.txt at rxn48573 — both immediately below the 2020 boundary at rxn48575. Unique_ModelSEED_Reaction_Aliases.txt and ..._Names.txt run all the way to rxn60859, so this is not missing source data: the annotation step simply stopped running. It was Scripts/Archived_Perl_Scripts/Compile_External_Pathways.pl, and it was never ported when the pipeline moved to Python.

A user asking "what pathway is this reaction in?" currently gets nothing for a fifth of the database.

What this PR does

Adds Scripts/Provenance/Build_Reaction_Pathways.py: joins ModelSEED → MetaCyc/KEGG alias → the committed pathway tables, walks the parent closure so class-level annotation is present, and writes both the alias file and the per-reaction pathways field.

What it deliberately does not do

It is additive, not a regeneration, and that is the whole design.

The committed pathway tables are a much smaller snapshot than the MetaCyc and KEGG distributions that produced the original alias file. Rebuilding every row from them reproduces only 104,203 of the 121,444 rows on disk and would silently drop 17,241 — every row of 601 reactions among them. Those distributions are licence-restricted and are not in this repository, so a faithful regeneration is not possible from a clean checkout and is not attempted.

So existing rows are left untouched and only reactions carrying no annotation are given one. 31,601 reactions remain unannotated, most of them the Rhea intake. Closing that needs the real pipeline port against the upstream sources; --check reports the remaining count so it stays visible rather than being quietly declared done:

$ python Build_Reaction_Pathways.py --check
alias rows on disk          121,444
  preserved                 121,444   (must equal the line above)
  added                      17,109
reactions annotated          21,734  ->  24,554
  newly annotated             2,820
still unannotated            31,601   of 56,012 reactions

EC numbers are not addressed. The same additive trick does not work: there is no committed EC table to join against, and Rhea's EC assignments need Data/rhea.rdf, which Scripts/Provenance/Rhea/Extract_Rhea_Reaction_Annotations.py expects and which is not in the repository. Left as a separate job rather than half-done.

Three things that had to be got right

Each was caught by checking the output, not by reasoning about it — worth recording since the next person to touch this will hit all three.

1. The flat files are not RFC-4180. The stoichiometry column carries unquoted " characters. A csv.writer round-trip quotes every such field and rewrote all 60 shards for a change touching a few hundred lines. The column is edited as text instead.

2. KEGG_pathways.tsv is a module table. Its ids are M-numbers and its parent column mixes rn pathway ids with free-text category names. Emitting the closure unfiltered produced entries like Metabolism (). The shipped annotation contains rn ids only, and so does this.

3. 683 of the 4,128 MetaCyc rows are ontology classes with an empty name — Degradation, Antibiotic-Biosynthesis, Fermentation. Emitting those directly would have written Antibiotic-Biosynthesis () 12,913 times. Display names already present in the alias file are reused instead, cutting it to 3,806 — the residue being classes whose name is recorded nowhere in the repository, which matches the 937 such rows the file already contained. No name is invented.

Verification

All checked before committing, on the full set:

• all 121,444 pre-existing alias rows preserved, none lost
• exactly 2,820 records changed, and the only field that differs on any of them is pathways
• no record that already had a value was overwritten
• every TSV keeps its line count (in-place column edit, no rows added or removed)
• the JSON shards differ only where a null became that record's list, so no shard is reformatted (60 shards, not 50 as first written)

Figure 1D in PR #299 does not depend on this PR — it joins the source table directly — so the two could have merged in either order.

🤖 Generated with Claude Code

Review pass (2026-09-23)

One defect: the 69 newly annotated records carrying both sources listed KEGG before MetaCyc; all 2,023 existing multi-source records are MetaCyc-first. Fixed — same 2,820 records, same rows, only element order on those 69. Re-verified: per-shard JSON/TSV change sets agree on all 50 shards; for spot-checked reactions (rxn50951, rxn58359, rxn60000) every direct source pathway is present in the emitted set with its ancestors.

Withdrawn (2026-09-28)

The maintainers' position is that pathways are not being included in the database — Biochemistry/REACTIONS.md already documents pathways as null — so this PR is closed. The branch is back on its standalone dev-based commits (a4039a2, 5da5c1e), which still merge cleanly onto dev; it had been rebuilt on the protonation-invariant work on 2026-09-25, and that work is now its own PR, #301. The branch stays on the fork should the position change. Nothing in #297–#299 depends on this: Figure 1D joins the MetaCyc source table directly.

One correction to the description above stands: a full rebuild from the committed tables reproduces 104,203 of the 121,444 alias rows, not 40,279, and would drop 17,241 (rxn00001 keeps all six of its own); the script's docstring on the branch still carries the old figure. The defect table is on the live basis; as first posted it divided live numerators by all-records denominators.